#identify modules

# Rscript step4.R --exprFile "../test_data/datExpr.csv" --network "../1_construct_network/network.RData" --output "../4_module_relationships"

library(optparse)


option_list = list(
  make_option("--exprFile", type="character", default=NULL,
              help="expression file path"),
  make_option("--network", type="character", default=NULL,
              help="network file path"),
  make_option("--output", type="character", default="out.txt",
              help="output path [default= %default]")
)

# 解析命令行参数
opt_parser = OptionParser(option_list=option_list, add_help_option=TRUE)
opts = parse_args(opt_parser)

datExpr = read.csv(opts$exprFile, stringsAsFactors = FALSE, row.names = 1);

library(WGCNA)

load(opts$network)
nGenes = length(net$goodGenes)

moduleColors = labels2colors(net$colors)


# Visualizing the network of eigengenes
png(
  filename = paste(opts$output, "module_relationships.png", sep="/"),
  width = 8,
  height = 16,
  units = "in",
  res = 300
)
MEs = moduleEigengenes(datExpr, moduleColors)$eigengenes
MET = orderMEs(MEs)
# Plot the relationships among the eigengenes and the trait
par(cex = 0.9)

plotEigengeneNetworks(
  MET,
  "",
  marDendro = c(0,4,1,2),
  marHeatmap = c(3,4,1,2),
  cex.lab = 0.8,
  xLabelsAngle = 90
)
dev.off()

# 计算模块间相关性
moduleCorrelation = cor(MEs, use = "p")
moduleCorPvalue = corPvalueStudent(moduleCorrelation, nrow(datExpr))

# 创建包含相关性和P值的矩阵
textMatrix = paste(signif(moduleCorrelation, 2), "\n(",
                  signif(moduleCorPvalue, 1), ")", sep = "")

# 计算FDR校正后的P值
moduleFDR = matrix(
    p.adjust(moduleCorPvalue, method = "BH"),
    nrow = nrow(moduleCorPvalue),
    dimnames = dimnames(moduleCorPvalue)
)

# 创建长格式的数据框，包含所有信息
result_df <- data.frame(
    Module1 = rep(rownames(moduleCorrelation), ncol(moduleCorrelation)),
    Module2 = rep(colnames(moduleCorrelation), each = nrow(moduleCorrelation)),
    Correlation = as.vector(moduleCorrelation),
    Pvalue = as.vector(moduleCorPvalue),
    FDR = as.vector(moduleFDR)
)

# 保存完整结果
write.csv(result_df,
          file = paste(opts$output, "module_relationships.csv", sep="/"),
          row.names = FALSE)
